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Abstract. We present results of a microlensing survey toward the Andromeda Galaxy (M31) carried out during four observing 
seasons at the Isaac Newton Telescope (INT). This survey is part of the larger microlensing survey toward M3 1 performed by 
the Microlensing Exploration of the Galaxy and Andromeda (MEGA) collaboration. Using a fully automated search algorithm, 
we identify 14 candidate microlensing events, three of which are reported here for the first time. Observations obtained at 
the May all telescope are combined with the INT data to produce composite lightcurves for these candidates. The results from 
the survey are compared with theoretical predictions for the number and distribution of events. These predictions are based 
on a Monte Carlo calculation of the detection efficiency and disk-bulge-halo models for M31. The models provide the full 
phase-space distribution functions (DFs) for the lens and source populations and are motivated by dynamical and observational 
considerations. They include differential extinction and span a wide range of parameter space characterised primarily by the 
mass-to-light ratios for the disk and bulge. For most models, the observed event rate is consistent with the rate predicted for 
self-lensing — a MACHO halo fraction of 30% or higher can be ruled at the 95% confidence level. The event distribution does 
show a large near-far asymmetry hinting at a halo contribution to the microlensing signal. Two candidate events are located 
at particularly large projected radii on the far side of the disk. These events are difficult to explain by self lensing and only 
somewhat easier to explain by MACHO lensing. A possibility is that one of these is due to a lens in a giant stellar stream. 

Key words. Gravitational lensing - M31: halo - Dark matter 



1. Introduction 

Compact objects that emit little or no radiation form a class of 
plausible candidates for the composition of dark matter halos. 
Examples include black holes, brown dwarfs, and stellar rem- 
nants such as white dwarfs and neutron stars. These objects. 
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collectively known as Massive Astrophysical Compact Halo 
Objects or MACHOs, can be detected indirectly through grav- 
itational microlensing wherein light from a background star is 
ampl ified by the spac e-time curvature associated with the ob- 
ject ( lPaczvnskJl986l) . 

The first microlensing s urveys wer e performed by the 
MACHO ('Alcock et al.'2000') and EROS /Lasserre et al."'200(^, 
Afonso et al. 2003) collaborations and probed the Milky Way 
halo by monitoring stars in the Large and Small Magellanic 
Clouds. While both collaborations detected microlensing 
events they reached different conclusions. The MACHO col- 
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laboration reported results that favour a MACHO halo fraction 
of 20%. On the other hand, the results from EROS are con- 
sistent with no MACHOs and imply an upper bound of 20% 
on the MACHO halo fraction. The two surveys are not in- 
consistent with each other since they probe different ranges in 
MACHO masses. They do leave open the question of whether 
MACHOs make up a substantial fraction of halo dark matter 
and illustrate an inherent difficulty with microlensing searches 
for MACHOs, namely that they must contend with a back- 
ground of self-lensing events (i.e., both lens and source stars 
in the Milky Way or Magellanic clouds), variable stars, and su- 
pernovae. The Magellanic Cloud surveys are also hampered by 
having only two lines of sight through the Milky Way halo. 

Microlensing surveys towards M3 1 hav e important a dvan- 
tages over the Magellanic Cloud surveys (Crotts 1992). The 
microlensing event rate for M31 is greatly enhanced by the 
high density of background stars and the availability of lines- 
of-sight through dense parts of the M31 halo. Furthermore, 
since lines of sight toward the far side of the disk pass through 
more of the halo than those toward the near side, the event 



a MACHO population should exhibit a 
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near-iar asymme tr 
iBaltz etal.l2003l. 

Unlike stars in the Magellanic Clouds, those in M31 are 
largely unresolved, a situation that presents a challenge for 
the surveys but one that can be overcome by a variety of 
techniques. To date microlensing events toward M31 have 
been reported by four different c ollaborations, VATT- Columbia 
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) and WeCAPP (Riffeser et al. 2003i). 



Recently, the POINT-AGAPE collaboration presented an 
analysis of data from three seasons of INT observations in 
which they concluded that "at least 20% of the halo mass 
in the direction of M31 m ust be in the form of MACHOs" 
dCalchi Novati et all l2005l) . Their analysis is significant be- 
cause it is the first for M3 1 to include a model for the detection 
efficiency. 

The MEGA collaboration is conducting a microlensing sur- 
vey in order to quantify the amount of MACHO dark matter in 
the M3 1 halo. Observations are carried out at a number of tele- 
scopes including the 2.5m Isaac Newton Telescope (INT) on La 
Palma, and, on Kitt Peak, the 1 .3m McGraw-Hill, 2.4m Hiltner, 
and 4m Mayall telescopes. The observations span more than 
4 seasons. The first three seasons of INT data were acquired 
jointly with the POINT- AGAPE collaboration though the data 
reduct ion and analys i s have been performed independently. 

In de Jong et al. ('2004') (hereafter Paper I) we presented 
14 candidate microlensing events from the first two seasons 
of INT data. The angular distribution of these events hinted 
at a near-far asymmetry albeit with low statistical significance. 
Recentlv ,An et al . (2004a) pointed out that the distribution of 
variable stars also shows a near-far asymmetry raising ques- 
tions about the feasibility of the M3 1 microlensing program. 
However, the asymmetry in the variable stars is likely caused 
by extinction which can be modelled. 

In this paper, we present our analysis of the 4-year INT 
data set. This extension of the data by two observing seasons 



compared to Paper I is a significant advance, but this data set 
is still only a subset of the MEGA survey. The forthcoming 
analysis of the complete data set will feature a further increase 
in time-sampling and baseline coverage and length. But there 
are more significant advances from Paper I. We improve upon 
the photometry and data reduction in order to reduce the num- 
ber of spurious variable-source detections. We fully automate 
the selection of microlensing events and model the detection 
efficiency through extensive Monte Carlo simulations. Armed 
with these efficiencies, we compare the sample of candidate 
microlensing events with theoretical predictions for the rate of 
events and their angular and timescale distributions. These pre- 
dictions are based on new self-consistent disk-bulge-halo mod- 
els ( Widrow & Dubinski 20o3) and a model for differential ex- 
tinction across the M3 1 disk. The models are motivated by pho- 
tometric and kinematic data for M31 as well as a theoretical 
understanding of galactic dynamics. 

Our analysis shows that the observed number of events can 
be explained by self-lensing due to stars in the disk and bulge 
of M3 1 , contrary to the findings of Calchi Novati et al. ( 20q1). 
Our results are consistent with a no MACHO hypothesis, al- 
though we cannot rule out a MACHO fraction of 30%. 

Data acquisition and reduction methods are discussed in 
Sect.|2] The construction of a catalogue of artificial microlens- 
ing events is described in Sect. |3l This catalogue provides the 
basis for a Monte Carlo simulation of the survey and is used, in 
Sect.|4] to set the selection criteria for microlensing events. Our 
candidate microlensing events are presented in Sect.|5] The ar- 
tificial event catalogue is then used in Sect.|6lto calculate the 
detection efficiency. Our extinction model is presented in Sect. 
In Sect.lSjthe theoretical models are described and the predic- 
tions for event rate and distribution are presented. A discussion 
of the results and our conclusions are presented in Sects. |5]and 

M 



2. Data acquisition and reduction 

Observations of M31 were carried out using the INT Wide 
Field Camera (WFC) and spread equally over the two fields 
of view shown in Fig. [2 The WFC field of view is approxi- 
mately 0.25n'' and consists of four 2048x4100 CCDs with a 
pixel scale of 0.333". The chosen fields cover a large part of 
the far side (SE) of the M31 disk and part of the near side. 
Observations span four observing seasons each lasting from 
August to January. Since the WFC is not always mounted on 
the INT, observations tend to cluster in blocks of two to three 
weeks with comparable-sized gaps during which there are no 
observations. 

Exposures during the first (1999/2000) observing season 
were taken in three filters, r', g' and i', which correspond 
closely to Sloan filters. For the remaining seasons (2000/01, 
2001/02, 2002/03), only the r' and i' filters were used. Nighffy 
exposure times for the first season were typically 10 minutes 
in duration but ranged from 5 to 30 minutes. For the remaining 
seasons the default exposure time was 10 minutes per field and 
filter. Standard data reduction procedures, including bias sub- 
traction, trimming and flatfielding were performed in IRAF. 
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Table 1. Overview of the number of epochs used for each field 
and fiher 
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57 
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27 


28 
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35 


32 


33 


30 


Total 


169 


169 


147 


138 



Fig. 1. the layout of the two INT Wide Field Camera (WFC) 
fields in M3 1 . A small part of the south field close to the bulge 
is not used since the image subtraction is not of high quality 
due to the high surface brightness. 

2.1. Astrometric registration and image subtraction 

We use Difference Image Photometry (DIP) 
(Tomaney & Crotts 1996) to detect variable objects in 
the highly crowded fields of M31. Individual images are 
subtracted from a high quality reference image to yield differ- 
ence images in which variable objects show up as residuals. 
Most operations are carried out with the IRAF package 
DIFIMPHOT. 

Images are transformed to a common astrometric reference 
frame. A high signal-to-noise (S/N) reference image is made by 
stacking high-quality images from the first season. Exposures 
from a given night are combined to produce a single "epoch" 
with Julian date taken to be the weighted average of the Julian 
dates of the individual exposures. 

Average point spread functions (PSFs) for each epoch and 
for the reference image are determined from bright unsaturated 
stars. A convolution kernel is calculated by dividing the Fourier 
transform of the PSF from an epoch by the PSF transform from 
the reference image. This kernel is used to degrade the image 
with better seeing (usually the reference image) before image 
subtraction is performed (Tomanev & Crotts 1996). 

Image subtraction does not work well in regions with very 
high surface brightness because of a lack of suitable, unsatu- 
rated stars. For this reason we exclude a small part of the south 
field located in a high-surface brightness region of the bulge 
(see Fig.Q. 

2.2. Variable source detection 

Variable sources show up in the difference images as residu- 
als which can be positive or negative depending on the flux of 
the source in a given epoch relative to the average flux of the 
source as measured in the reference image. However, difference 
images tend to be dominated by shot noise. The task at hand is 
to differentiate true variable sources from residuals that are due 

to noise. 

The program SExtractor dBertin & Arnoutsll 19961) is used 
to detect "significant residuals" in r' epochs, defined as groups 
of 4 or more connected pixels that are all at least 3cr above 



or below the background. Residuals from different epochs are 
cross-correlated and those that appear in two or more consec- 
utive epochs are catalogued as variable sources. (Because of 
fringing, the i' difference images are of poorer quality than the 
r' ones and we therefore use r' data to make the initial identifi- 
cation.) 

2.3. Lightcurves and Epoch quality 

The difference images for a number of epochs are discarded 
for a variety of reasons. Epochs with poor seeing do not give 
clean difference images. We require better than 2" seeing and 
discard 7 epochs and parts of 12 epochs where this condition is 
not met. PSF-determination fails if an image is over-exposed. 
We discard 7 epochs and parts of another 7 epochs for this rea- 
son. Finally 2 epochs from the second and third seasons are 
discarded because of guiding eiTors. 

Lightcurves for the variable sources are obtained by per- 
forming PSF-fitting photometry on the residuals in the differ- 
ence images. For every pixel the Poisson-noise is evaluated as 
well as the fractional flux eiTor due to photometric inaccura- 
cies in the matching and subtraction steps for the difference 
image in question. Fluxes and their error bars are derived by 
optimal weighting of the individual pixel values. Lightcurves 
are also produced at positions where no variability is identi- 
fied and fit to a flat line. These lightcurves serve as a check 
on the contribution to the flux error bars derived from the pho- 
tometric accuracy of each difference image. For each epoch, 
we examine by eye the distribution of the deviations from the 
flat-line fits normalised by the photometric eiTor bar Epochs 
where this distribution shows broad non-Gaussian wings are 
discarded since wings in the distribution are likely caused by 
guiding errors or highly variable seeing between individual ex- 
posures. For epochs where this normalised error distribution is 
approximately Gaussian but with a dispersion greater than one, 
the error bars are renormalised. 

Approximately 19% of the 209 r' epochs and 22% of the 
183 i' epochs are discarded. The number of epochs that remain 
for each season, filter, and field are tabulated in Table[2 Though 
variable objects are detected in r', lightcurves are constructed 
in both r' and i'. In total, 105,447 variable source lightcurves 
are generated. 

3. Artificial microlensing events 

This section describes the construction of a catalogue of ar- 
tificial microlensing lightcurves which forms the basis of our 
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Monte Carlo simulations. We add artificial events to the differ- 
ence images and generate lightcurves in the same manner as is 
done with the actual data. The details of this procedure follow 
a review of microlensing basics and terminology. 

3.1. Microlensing lightcurves 

The lightcurve for a single-l ens microlensing event is described 
by the time-dependent flux (|Paczvriskill986l) : 



where 



F(t) = Fo 



+2 



u Vm^ -1-4 



= FoA(f) 



(1) 



where F{) is the unlensed source flux and A is the amplification. 
u - u(t) is the projected separation of the lens and the source 
in units of the Einstein radius. 
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where m is the lens mass and the D's are the distances between 
observer, lens and source. If the motions of lens, source, and 
observer are uniform for the duration of the lensing event we 
can write 



(3) 



where fi is the impact parameter in units of R^, that is, the mim- 
imum value attained by u. fmax is the time of maximum ampli- 
fication and fE is the Einstein time, defined as the time it takes 
the source to cross the Einstein radius. 

In classical microlensing the measured lightcurves contain 
contributions from unlensed sources. Blending, as this effect is 
known, changes the shape of the lightcurve and can also spoil 
the achromaticity implicit in equation^ In our survey, we mea- 
sure flux differences that are created by subtracting a reference 
image. Since the flux from unlensed sources is subtracted from 
an image to form the difference image, blending is not a prob- 
lem unless the unlensed sources are variable. Blending by vari- 
able sources does introduce variations in the baseline flux and 
adversely affects the fit. 

For a difference image the microlensing lightcurve takes the 
form 



mt) = Fit) - F,ef = AFbi + Fo(A(0 - 1) 



(4) 



where fref is the reference image flux and Afbi = - ^ref- 
Thus, if in the reference image the source is not lensed, F^^i - 
Fo and therefore AF\,i = 0. Only if the source is amplified in 
the reference image will AFbi be non-zero and negative. 

For unresolved sources, a situation known as pixel lensing 
(and the one most applicable to stars in M31), those microlens- 
ing events that can be detected typically have high amplifica- 
tion. In the high amplification limit, Ie and /? are highly degen- 
erate (lOould 1996; Baltz & Silk 2000) and difficult to extract 
from the lightcurve. It is therefore advantageous to parame- 
terise the event duration in terms of the half-maximum width 
of the peak. 



and 

^Jx{x + 4) 



(6) 



(7) 



fpWHM = fEW(/3) 



(5) 



dGondolol 19991) . w(J3) has the limiting forms wifi -^l) - fi^ 
andw08» \) ^ - 1^1^. 

3.2. Simulation parameters 

The parameters that characterise microlensing events fall into 
two categories: "microlensing parameters" such as /3, fmax, 
and tE, and parameters that describe the source such as its 
brightness Fq i-, its r'-i' colour C, and its position. We survey 
many lines-of-sight across the face of M31. Furthermore, 
all types of stars can serve as a source for microlensing. 
Therefore, our artificial event catalogue must span a rather 
large parameter space. This parameter space is summarised in 
Tableland motivated by the following arguments: 

• Peak times and baseline fluxes 

We demand that the portion of the lightcurve near peak 
amplitude is well-sampled and therefore restrict f^ax to one 
of the four INT observing seasons. The reference images are 
constructed from exposures obtained during the first season. If 
a microlensing event occurs during the first season and if the 
source is amplified in one or more exposures during this sea- 
son, the baseline in the difference image will be below the true 
baseline. For an actual event in season one, this off-set is ab- 
sorbed in one of the fit parameters for the lightcurve. For artifi- 
cial events, the baseline is coiTected by hand. 

• Event durations 

Limits on the duration of detectable events follow naturally 
from the setup of the survey and the requirement that events 
are sampled through their peaks. Since the INT exposures are 
combined nightly, events with fpwHM < 1 day are practically 
undetectable except for very high amplifications. On the other 
hand, events with fpwHM approaching the six-month length of 
the observing season are also difficult to detect with the selec- 
tion probability decreasing linearly with frwHivi- Because gaps 
in the time coverage of our survey will affect our sensitivity to 
short events more strongly than to long events, sampling should 
be denser at shorter timescales. To limit computing time and 
ensure statistically significant results spread over a wide range 
of event durations, we simulate events at six discrete values of 
tpwHM' 1, 3, 5, 10, 20 and 50 days. 

• Source fluxes and colours 

Faint stars are more abundant than bright ones. On the other 
hand, microlensing events are more difficult to detect when the 
source is a faint star. The competition between these two ef- 
fects means that there is a specific range of the source lumi- 
nosity function that is responsible for most of the detectable 
microlensing events. 



Jelte T.A. de Jong et al.: MACHOs in M31? Absence of evidence but not evidence of absence 



5 



-5 



-2 




10-3 0.01 0.1 1 10 

r' flux (ADU/s) 



100 



Fig. 2. The solid line in this figure shows the R-band luminos- 
ity function from Mamon & Soneira (1982). Multiplying this 
function with the square of the maximum impact parameter 
y6max needed to detect a microlensing event gives the dashed 
line. The line shown is for a detection threshold of 1 ADUs ' 
in r'. The upper horizontal axis shows absolute R-band magni- 
tude, the lower axis the corresponding r' flux. 



The maximum flux difference during a microlensing event 



IS 



AFrr 



- 1 



(8) 



where we are ignoring the AF\,i term in equation |3 Let AFdet 
be the detection threshold for AF^ax- A lower bound on AFniax 
implies an upper bound on /3 which, through equation |8j is a 
function of the ratio Fo/Fdei- j6u = Pu (Fo/Fdi,,). The probabil- 
ity that a given source is amplified to a detectable level scales 
as ySy. In Fig. El we show both th e R-ba nd luminosity func- 
tion, A^,, from lMamon & Soneiral ( Il982l) and the product of 
this luminosity function with fi^ assuming a detection thresh- 
old of Fdet - 1 ADU s The latter provides a qualitative pic- 
ture of the distribution of detectable microlensing events. This 
distribution peaks at an absolute R-band magnitude of approx- 
imately indicating that most of the sources for detectable mi- 
crolensing events are Red Giant Branch (RGB) stars. 

Since there is no point in simulating events we cannot de- 
tect we let the impact parameter /3 vary randomly between 
and ySu- Table|2lsummarises the fluxes and values for /3u used in 
the simulations. 

For the artificial event catalogue, we use source stars with 
a r' fluxes at several discrete values between 0.01 and 10 
ADUs '. Typically the r'-i' colours of RGB stars range be- 
tween C = 0.5 and 2.0. We assume C = 0.75 for our artificial 
events. As a check of the dependence of the detection efficiency 
with colour, we also simulate events with C = 1.25. 
• Position in M31 



Table 2. Fluxes and maximum impact parameters probed in the 
simulations of microlensing events. 
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0.01 


29.5 


0.011 


28.75 


0.01 


0.1 


27.0 


0.11 


26.25 


0.09 


0.5 


25.2 


0.55 


24.45 


0.35 


1.0 


24.5 


1.11 


23.75 


0.56 


10.0 


22.0 


11.1 


21.25 


1.67 



Lightcurve quality and detection efficiency vary with po- 
sition in M31 for several reasons. The photometric sensitivity 
and therefore the detection efficiency depend on the amount of 
background light from M31 and are lowest in the the bright 
central areas of the bulge. Diff'erence images from these areas 
are also highly crowded with variable-star residuals which 
influence the photometry and add noise to the microlensing 
lightcurves. To account for the position-dependence of the 
detection efficiency, artificial events are generated across the 
INT fields. To be precise, the artificial event catalogue is 
constructed in a series of runs. For each run, artificial events 
are placed on a regular grid with spacing of a 45 pixels (^15 ") 
so that there are 3916 artificial events per chip. The grid is 
shifted randomly between runs by a maximum of 10 pixels. 

To summarise, artificial events are characterised by the pa- 
rameters fpwHM, ^0, C, tma^,/3, and their angular position. These 
events are added as residuals to the difference images using the 
PSF in the subregion of the event. The residuals also include 
photon noise. The new difference images are analysed as in 
Sect.lJland lightcurves are built for all artificial events detected 
as variable objects. 

4. Microlensing event selection 

The vast majority of variable sources in our data set are vari- 
able stars. In this section we describe an automated algo- 
rithm that selects candidate microlensing lightcurves from this 
rather formidable background. Our selection criteria pick out 
lightcurves that have a flat baseline and a single peak with the 
"correct" shape. The criteria take the form of conditions on the 

statistic that measures the goodness-of-fit of an observed 
lightcurve to equation^] The fit involves seven free parameters: 
fmax,y6, fE , F()j, F(),i, AFbi i-, and AFbi.i. To increase computation 
speed we first obtain rough estimates for fmax and Ie from the r' 
lightcurve and then perform the full 7-parameter fit using both 
r' and i' lightcurves. 

Gravitational lensing is achromatic and therefore the ob- 
served colour of a star undergoing microlensing remains con- 
stant in contrast with the colour of certain variables. While we 
do not impose an explicit achromaticity condition, changes in 
the colour of a variable source show up as a poor simultane- 
ous r' and i' fit. Because many red variable stars vary little in 
colour, as defined by measurable differences in flux ratios, the 
lightcurve shape and baseline flatness are better suited for dis- 
tinguishing microlensing events from long period variable stars 
(LPVs) than a condition on achromaticity. 
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Fig. 3. Scatter plots of t^x^ vs. for simulated events with fFWHM=50 days (a), 10 days (b), 1 day (c), and for the actual data 
for 1 CCD. The solid lines correspond to equations 9 and 10. 



Lightcurves must contain enough information to fit ade- 
quately both the peak of the microlensing event and the base- 
line. We therefore impose the following conditions: (1) The r' 
and i' lightcurves must contain at least 100 data points. (2) The 
peak must be sampled by several points well-above the base- 
line. (3) The upper half of the peak, as defined in the difference- 
image lightcurve, must lie completely within a well-sampled 
observing period. The second condition can be made more pre- 
cise. We allow for one of the following two possibilities: (a) 4 
or more data points in the r'-lightcurve are 3cr above the base- 
line or (b) 2 or more points in r' and 1 or more points in i' are 
3cr above the baseline. (The r' data is weighted more heavily 
than the i' data because it is generally of higher quality and be- 
cause i' was not sampled as well during the first season.) The 
third condition insures that we sample both rising and falling 
sides of the peak. We note that there are periods during the last 
two seasons where we do not have data due to bad weather. 
The periods we use ai-e the following: 01/08/1999-13/12/1999, 
04/08/2000-23/01/2001, 13/08/2001-16/10/2001, 01/08/2002- 
10/10/2002, and 23/12/2002-31/12/2002. 

The selection of candidate microlensing events is based on 
the ;t'^-statistic for the fit of the observed lightcurve to equa- 
tionl^as well as t^x^ = x\^i — '*here;^fg^^ is the;t'^-statistic 
for the fit of the observed lightcurve to a flat line. Our^^-cuts 
are motivated by simulations of artificial microlensing events. 



In Fig. |3] we show the distribution of artificial events with 
fpWHM = 50, 10, and 1 days (panels a, b, and c respectively) 
and for all variable sources in one of the CCDs (panel d). In 
Fig. 15 we show the variable sources from all CCDs that satisfy 
conditions 1-3. The plots are presented in terms of x~ IN and 
Ax^/N where is the number of data points in an event. We 
choose the following cuts: 

AX^ > 1.5N (9) 
and 

X^ <(N-l)f{Ax'-) + 3(2(N-7)y'^ (10) 

where / ( A^^) = Ax^ / 1 00 + 1 . The first criterion is meant to fil- 
ter out peaks due to noise or variable stars. The second criterion 
corresponds to a 3cr-cut in x^ for low signal-to-noise events. 
The x^ threshold increases with increasing Ax^ ■ Panels a-c of 
Fig.Elshow a trend where x^ increases systematically with Ax^- 
This effect is due to the photometry routine in DIFIMPHOT 
which underestimates the error in flux measurements for high 
flux values. The function / is meant to compensate for this ef- 
fect. 

The selection criteria appear as lines in Figs.|3and|3 (To 
draw these lines, we take A^ - 309 though in practice N is 
different for individual lightcurves.) 
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Fig. 4. Ax^/N versus /N for variable sources that satisfy se- 
lection criteria (1), (2) and (3) for peak and lightcurve sam- 
pling. The solid line indicates criteria (4) and (5) for peak sig- 
nificance and goodness of fit. Criterium (5) depends on the 
number of points in the lightcurves, and the line drawn here 
is for A^=309, the typical number of available data points per 
source. Two candidate events with higher Ax^/N are indicated 
with arrows, labelled with their Ax^/N value. 

5. Candidate events 

Of the 105 477 variable sources 28 667 satisfy conditions 1-3. 
Of these, 14 meet the criteria set by equations 1^ and [TUI The 
positions of 12 of these events in thex^/N - Ax^/N plane are 
shown in Fig.|3 

5.1. Sample description 

In Table |3] we summarise the properties and fit parameters of 
the 14 candidate microlensing events. The first column gives 
the assigned names of the events using the nomenclature from 
Paper 1. The numbering reflects the fact that candidates 4, 5, 

6, and 12 from Paper 1 are evidently variable stars since they 
peaked a second time in the fourth season. The other 10 events 
from Paper 1 are "rediscovered" in the cuiTent more robust 
analysis. Four additional candidates, events 15, 16, 17, and 
18, are presented. Even t 16 is the same as PA-99-N1 from 
IPaulin-H enriksson et al.l ( l2003l) and was not selected in our 
previous analysis because the baseline was too noisy due to a 
nearby bright variable star. It now passes our selection criteria 
thanks to the smaller aperture used for the photometry (see dis- 
cussion below). The three other events all peaked in the fourth 
observing season and are reported here for the first time. 

The coordinates of the events are given in columns 2 and 3 
of Table |3j their positions within the INT fields are shown in 
Fig- HI The fit parameters, and Ax^ are given in the remain- 
ing columns. In AppendixIXlwe show the r' and i' lightcurves, 
thumbnails from the difference images for a number of epochs. 



fig05.jpg 



Fig. 5. The locations of the 14 microlensing events within the 
INT fields are shown here with the d ots. Events 7 and 16 cor- 
respon d with events N2 and Nl from lPaulin-Henriksson et all 
(I2OO3I) . Their event S3 is indicated with a cross and lies in the 
high surface brightness region that we exclude from our anal- 
ysis. Also ma rked with a cross (B 1 ) is the position of level 1 
candidate 1 of lBelokurov et al.l ( l2005h . 

and a comparison of Ar' and Ai' for points near the peak. The 
latter provides an indication of the achromaticity of the event. 
The lightcurves include data points from observations at the 4m 
Mayall telescope on Kitt Peak (KP4m) though the fits use only 
INT data. 

We have already seen that variable stars can mimic mi- 
crolensing events. Blending of variable stars is also a problem 
since it leads to noisy baselines. This problem was rather se- 
vere in Paper 1 causing us to miss event PA-99-N1 found by the 
POINT- AG APE collaboration ( Paulin-Henriksson et al. 20031) . 
In an eff'ort to reduce the effects of blending by variable stars, 
we use a smaller aperture when fitting the PSF to residuals in 
the difference images. Nevertheless, some variable star blend- 
ing is unavoidable, especially in the crowded regions close to 
the centre of M3 1 . Event 3 provides an example of this effect. 
A faint positive residual is visible in the 1997 KP4m differ- 
ence image as shown in Fig.|6| The residual is located one pixel 
(0.21") from the event and is likely due to a variable star. It cor- 
responds to the data point in the lightcurve ~1000 days before 
the event and well-above the baseline (see Fig. lA.3t . The KP4m 
data point from 2004 is also above the baseline but in this and 
other difference images, no residual is visible. The implication 
is that variable stars can influence the photometry even when 
they are too faint to be detected directly from the difference 
images. 

Good simultaneous fits are obtained in both r' and i' for all 
candidate events. Event 7 has a high;if^/A^ of 1.98, but since 
Ax^/N is very high, the event easily satisfies our selection cri- 
teria. In high S/N events, secondary effects from parallax or 
close caustic approaches can cause measurable deviations from 
the standard microlensing fit. In addition, as discussed above, 
we te nd to underestima te the photometric errors at high flux 
levels. lAnet alJ (12004b') studied this event in detail and found 
that the deviations from the standard microlensing shape of the 
POINT- AGAPE lightcurve are best explained by a binary lens. 
The somewhat high;^'^ for events 10 and 15 are probably be- 
cause they are located in regions of high surface brightness. 



8 



Jelte T.A. de Jong et al.: MACHOs in M31? Absence of evidence but not evidence of absence 



Table 3. Coordinates, highest measured difference flux, and some fit parameters for the 14 candidate microlensing events. The 
peak time f^ax is in days after August 1st 1999 (JD 2451393). 



Candidate 


RA 


DEC 


Ar' 






fpWHM 


X IN 




^0,r 


r'-i' 


event 


/ T'^ r\t\f\\ 

(J2000) 


(J 2000) 


(mag) 


(days) 


(days) 






(ADU s"') 


(mag) 


MEGA-ML 1 


0:43:10.54 


41:17:47.8 


21.8±0.4 


60.1 


±0.1 


5.4 


±7.0 


1.12 


1.91 


0.1±0.3 


0.6 


MEGA-ML 2 


0:43:11.95 


41:17:43.6 


21.51±0.06 


34.0 


±0.1 


4.2 


±0.7 


1.06 


2.48 


3.4±1.7 


0.3 


MEGA-ML 3 


0:43:15.76 


41:20:52.2 


21.6±0.1 


420.03 4 


0.03 


2.3 


±2.9 


1.14 


2.11 


0.08±0.21 


0.4 


IVliiOA-iVIL / 




41 .Zo. 44.0 


Ly.j I +U.UZ 


71.8 


±0.1 


17.8 


±0.4 


1 QQ 

V.yo 




0.o±U.4 


1 c 
L.J 


MEGA-ML 8 


0:43:24.53 


41:37:50.4 


22.3±0.2 


63.3 


±0.3 


27.5 


± 1.2 


0.82 


3.03 


20.4±22.9 


0.6 


MEGA-ML 9 


0:44:46.80 


41:41:06.7 


21.97±0.08 


391.9 


±0.1 


2.3 


±0.4 


1.02 


2.49 


0.9±0.4 


0.2 


MEGA-ML 10 


0:43:54.87 


41:10:33.3 


22.2±0.1 


75.9 


±0.4 


44.7 


±5.6 


1.28 


5.88 


1.4±0.5 


1.1 


MEGA-ML 1 1 


0:42:29.90 


40:53:45.6 


20.72±0.03 


488.43 4 


0.04 


2.3 


±0.3 


1.03 


13.27 


1.5±0.4 


0.2 


MEGA-ML 13 


0:43:02.49 


40:45:09.2 


23.3±0.1 


41.0 


±0.3 


26.8 


± 1.5 


0.75 


1.68 


9.2±10.8 


0.8 


MEGA-ML 14 


0:43:42.53 


40:42:33.9 


22.5±0.1 


455.9 


±0.1 


25.4 


±0.4 


1.11 


3.74 


146±182 


0.4 


MEGA-ML 15 


0:43:09.28 


41:20:53.4 


21.63+0.08 


1145.5 


±0.1 


16.1 


± 1.1 


1.23 


4.41 


7.0±2.2 


0.5 


MEGA-ML 16 


0:42:51.22 


41:23:55.3 


21.16±0.06 


13.38 4 


0.02 


1.4 


±0.1 


0.93 


2.81 


2.6±0.7 




MEGA-ML 17 


0:41:55.60 


40:56:20.0 


22.2±0.1 


1160.7 


±0.2 


10.1 


± 2.6 


0.79 


2.02 


0.5±0.3 


0.4 


MEGA-ML 18 


0:43:17.27 


41:02:13.7 


22.7±0.1 


1143.9 


±0.4 


33.4 


±2.3 


1.13 


1.83 


13.7±16.3 


0.5 




Fig. 6. Detail of two KP4m difibrence images centred on the 
position of event 3. Left: October 27th 1997, almost 3 years 
before the event peaks, a very faint residual is seen centred just 
1 pixel (0.21") away from the event. Right: September 26th 
2000, during the peak of the event that is displaced from the 
position of the faint variable. 



All of the candidate events are consistent with achromatic- 
ity, though for events with low S/N, it is difficult to draw firm 
conclusions directly from the lightcurves or Ar' vs. A/' plots. 
The values for Fo.r and C for the events give some indication 
of the properties of the source stars. The unlensed fluxes are 
consistent with the expected range of 0. 1 - 10 ADU s ' and the 
colours for most of the events are typical of RGB stars. Note 
however that for many of the events, the uncertainties for Fq, 
/3, and fpwHivi are quite large. These uncertainties reflect degen- 
eracies among the lightcurve fit parameters. 

The number of candidate events varies considerably from 
season to season. We find 7 events in the first season, 4 in the 
second season, none in the third season and 3 in the fourth sea- 
son. The paucity of events during the third and fourth seasons 
is not surprising given that we have fewer epochs for those 
seasons (see table In particular, the gaps in time coverage 
during those seasons conspired against the detection of short 
duration events. 
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Fig. 7. Our photometry for microlensing event candidate 1 from 
iRelokiirovetalllEooi. 



5.2. Comparison with other surveys 

The POINT- AGAPE collabor ation published several analyses 
of the INT observations. In IPaulin-Henriksson et alJ (l2003h 
they presented four convincing microlensing events from the 
first two observing seasons using stringent selection criteria. In 
particular, they restricted their search to events with high S/N 
and fpwHM < 25 days. They argued that one of these events 
(PA-00-S3) is probably due to a stellar lens in the M31 bulge. 
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log(Fsource (ADU/s)) 

Fig. 8. Relative probability of detecting a microlensing event 
of a source star with a certain intrinsic flux. This probability 
is the product of the number of available stars (taken from the 
luminosity function), the square of the maximum impact pa- 
rameter for which an event can be detected, and the detection 
efficiency for each source population, averaged over all fpwHM- 

This event lies in the region of the bulge excluded from our 
analysis (see Fig.^J. The other three events, PA-99-N1, PA-99- 
N2, and PA-00-S4, correspond respectively to our events 15, 7, 
and 1 1 . Evidently, the remaining eight events from our analysis 
of the first two INT seasons did not satisfy their rather severe 
select ion criteria. 

In lBelokurov et alJ ll2005h . flie POINT- AGAPE coUabora- 
tion analysed data from the first three INT observing seasons 
without any restrictions on the event duration. Using differ- 
ent selection criteria from their previous analyses, they found 
three high quality candidates. Two of these events were already 
known (PA-00-S4 or MEGA-ML-11 and PA-00-S3). The one 
new event is present in our survey but does not pass our se- 
lection criteria because of a high ■ The lightcurve for this 
event, along with our best-fit model, is shown Fig. The 
model does not do a good job of reproducing the observed 
lightcurve behaviour In particular, the observed lightcurve ap- 
pears to be asymmetric about the peak time f^ax- The observed 
r'-lightcurve is systematically below the model 15-20 days 
prior to f^ax- Both r' and i' lightcurves are above the model 
10-15 days after f^ax- Since there are no data available on the 
rising part of the peak, f^ax is poorly constrained and may in 
fact be less than the 770 days used in the fit. The shape of our 
x' ligh tcurve is similar to the one presented in .Belokurov et al. | 
( I2OO5I) (NB. They removed one epoch close to the peak cen- 
tre that is present in our lightcurve.) In i' the peak shapes are 
somewhat different. 

Peak asymmetries can be caused by secondary effects such 
as parallax. In our opinion, a more likely explanation for this 
case is that the event is a nova-like eruptive variable. Granted, 
the event appears to be achromatic. But classical novae can be 
achro matic on the declining part of the lightcurve (see, for ex- 
ample, jDamle^raL, (.2004.) ). precisely where there is data. If 



this is a classical nova, it would be a very fast one, with a de- 
cli ne rate corresponding to ~0 .6 mag per day. 

ICalchi Novati et alT (l2005) found six candidate microlens- 
ing events in an analysis of the three-year INT data 
set. Of these events, four a re the same as reported by 
IPauUn-Henriksson et alJ ( l2003l) and two are new events: PA- 
00-N6 and PA-99-S7. The latter of these is located in the bright 
part of the southern field excluded in our analysis (Fig. Q. 
Candidate event PA-00-N6 is present in our data, but was only 
detected in one epoch in our automatic SExtractor residual de- 
tection step and t herefore did not make it in to the catalogue of 
variable sources. ICalchi Novati et alJ (12005') do not detect our 
events 1, 2, 3, 8, 9, 10, 13, and 14, which all peak in the first 
two observing seasons. Evidently, these events do not satisfy 
their S/N constraints. 

6. Detection efficiency 

We determine the detection efficiency for microlensing events 
by applying the selection criteria from Sect.|31to the catalogue 
of artificial events from Sect.|3l As discussed above, simulated 
lightcurves are generated by adding artificial events to the dif- 
ference images and then passing the images through the pho- 
tometry analysis routine designed for the actual data. Those 
lightcurves that satisfy the selection criteria for microlensing 
form a catalogue of simulated detectable microlensing events. 
The detection efficiency is the ratio of the number of these 
events to the original number of artificial events. 

We first check that our artificial event catalogue includes 
the portion of the source luminosity function responsible for 
most of the detectable events. The function N»ii\ in Fig.|2]is 
meant to give a qualitative picture of the detectability of mi- 
crolensing as a function of source luminosity. Here we con- 
sider the function Pdet = Ntfile where e is detection efficiency 
as a function of ^ integrated over fi, fpwHM and position, f det 
gives the relative probability for detection of a microlensing 
event as a function of the source luminosity. As shown in Fig. 
|SJ the range 0.01 to 10 ADU s~' adequately covers the peak of 
this probability distribution. 

Our goal is to represent the detection efficiency in terms of 
a simple portable function of a few key parameters. We adopt a 
strategy whereby the detection efficiency is modelled as func- 
tions of fpwHM and AF^ax for individual subregions of the two 
fields. The parameters p and fmax are "integrated out" and C is 
fixed to the value 0.75. This strategy is motivated by the fol- 
lowing considerations. 

In Fig.|9]we plot the detection efficiencies as a function of 
P for four different values of fpwHM- In each of the panels, the 
efficiencies are integrated over position within a single chip of 
the INT fields. The top (bottom) panels are for the south-east 
chip of the north (south) field. The right (left) panels are for 
bright (faint) source stars. The general trend is for the detec- 
tion efficiency to increase with increasing fpwHM and decreas- 
ing p. This trend is expected since longer duration events are 
more likely to be observed near the peak and smaller values of 
P imply larger amplification factors. For Fir) - 10 ADU s ', 
fpWHM ^ 10 days and small p < 0.7, the detection efficiencies 
decrease with decreasing p. The decrease is more severe for 
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Fig. 9. Detection efficiencies as function of impact parameter p for different values of fpwHM (50, 10, 3 and 1 days). The two 
upper panels show the fraction of simulated events that pass the microlensing selection criteria for 2 source fluxes, 10 and 0.01 
ADU s in the south-east chip of the north field. The lower panels show the same for the south-east chip of the south field. 



the fpwHM = 50 day events where the detection efficiency ac- 
tually drops below that for the fpwHM - 10 day events. The 
problem may be that we underestimate the photometric error 
at high fluxes therefore causing to be systematically high. 
Moreover, 50 days is a substantial fraction of the observing 
season and therefore some long duration events may not meet 
the requirement that the peak be entirely within a single season. 

Since the shape of the microlensing lightcurve does not 
depend strongly on (3 we expect no significant dependence of 
the detection efficiency on the intrinsic source brightness. This 
point is illustrated in Fig. 1 101 where we plot the detection effi- 
ciencies as a function of l/AF^ax for events with fpwHM = 50 
days. We integrate the efficiencies over positions within single 
CCDs and show the results for four of the eight CCDs in our 
fields. The curves vary by at most 30% over three orders of 
magnitude in F(r). The implication is that an explicit F(r) de- 
pendence in the detection efficiency will not change the results 
significantly. 

We next test whether the detection efficiency depends on 
the colour C of the source. In addition to the main artificial 
event catalogue, we generate artificial events with C - 1 .25 and 
r' unchanged for a part of ffie north field. Fig.^Jcompares the 
detection efficiencies for the two colours and shows that there 
is no significant difference, except for the very highest signal to 



noise events. The discrepancy at high S/N reflects the problem 
discussed above with our estimates of the photometric errors at 
high flux. This problem is worse for redder sources which have 
a higher i'-band flux. 

Motivated by the shapes of the curves in Fig.[n] we choose 
a Gaussian in l/AFmax where the position of the peak depends 
on fpwHM- The explicit functional form is taken to be: 

e = ci(l-fpwHM/112)e-^('/^^"«-^-'^' (11) 
where 

C3 = di ■In(fpwHM) + di ■ (12) 

The factor multiplying the Gaussian takes into account the 
sharp decrease in detection efficiency for events with duration 
comparable to or longer ffian the observing season. The param- 
eters ci,C2,d\ and d2 are determined by fitting simultaneously 
the detection efficiencies for all values of fpwHM to equation 
[TTI Fig. shows an example of ffiese fitting formulae to the 
detection efficiencies. 

Fig.^Jillustrates the dependence of the detection efficien- 
cies on location in the INT fields. This dependence is due 
mainly to variations in galaxy surface brightness but also to the 
presence of bad pixels and saturated- star defects. As discussed 
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Fig. 10. Detection efficiencies as a function of l/AF^^j^ for fpwHM - 50 days and ^ =10 ADUs ' (solid line), 1 ADUs ' 
(dotted line), 0.1 ADU s"' (long-dashed), and 0.01 ADU s ' (short-dashed line). In general the lines overlap within the errors. 




0.5 



Fig. 11. Colour dependence of the detection efficiency. For 
fpWHM 's of 1, 3, 10 and 50 days the detection efficiencies are 
shown for the 2 different source colours simulated. The colour 
has no noticeable effect, except for the highest signal-to-noise 
events. 



above, we account for the spatial dependence by fitting the de- 
tection efficiency separately for subregions of the fields. To be 
precise, we divide each chip into 32 subregions, ~3'x3' in size. 
For each of these regions we average 14640 simulated events 
(2440 per choice of fpwHM)- 



7. Extinction 

Microlensing surveys such as MEGA and POINT-AGAPE are 
motivated, to a large extent, by the argument that a MACHO 
population in M3 1 would induce a near-far asymmetry in the 
microlensing event distribution. In the absence of either extinc- 
tion or significant intrinsic asymmetries in the galaxy, the dis- 
tribution of self-lensing events and variable stars masquerading 
as microlensing events would be near-far symmetric. The de- 
tection of a near-far asymmetry would then provide compelling 
evidence in favour of a significant MACHO population. 

Recently, lAn et alJ ll2004al) found a near-far asymmetry in 
the distribution of variables which they attribute to differential 
extinction across the M31 disk. That differential extinction is 
significant is also witnessed by several dust features including 
two prominent dust lanes on the near side of the disk. 

We construct a simple model for differential extinction in 
M3 1 and test it to against the distribution of LPVs. In the next 
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Fig. 12. Detection efficiencies as a function of 1 / AF^^^ for 
different values of fpwHM- The symbols give the results of the 
Monte Carlo calculation for one chip. The lines correspond to 
the fitting formula, eauationFTTI 



section, we incorporate this extinction model into our calcula- 
tions for the t heoretical event rate. 

Following IWalterbos & KennicuttI lll98 8l) we assume that 
the dust is located in a thin layer in the mid-plane of the disk. 
Along a given line-of-sight, only light from behind the dust 
layer is absorbed. Because of the galaxy's high inclination, the 
fraction of stars located behind the dust layer is higher for lines- 
of-sight on the near side of the disk than for those on the far 
side, as illustrated in Fig.^] Therefore, even if the distribution 
of dust is intrinsically symmetric, extinction will have a greater 
effect on the near side of the disk. 

Based on these assumptions the observed intensity along a 
particular line-of-sight is 



I. 



obs 



-^front + i'backe 



(13) 



where Jfiont (J^hack) is the intensity of light originating from in 
front of (behind) the dust layer and t is the optical depth. This 
equation can be rewritten in terms of the total intrinsic intensity, 
Jintr, and the fraction x of light that originates from in front of 
the dust layer: 



J, 



obs 



JcJintr + (1 - X)linu-e 



(14) 



The three unknowns in this equation, lina, x, and e depend 
on wavelength. Rewriting equation^Jfor the B-band we have 



Jobs(g)/Jintr(g) - XB 
1 - Xb 



(15) 



As a first approximation we assume that Jobs(^) = -^intr(^) so 
that 



i,M/(Cbi ■ -robs(/)) 

1 -X 



(16) 




Near 



Far 



Fig. 13. Schematic representation of the line-of-sight through 
the M3 1 galaxy from an observer on earth. Because of the high 
inclination of M31, most of the light observed on the near side 
of the disk is coming from behind the dust lanes. 

Table 4. Disk and bulge parameters used to derive x, the frac- 
tion of light originating in front of the mid-plane of M3 1 : the 
scale length and scale height, hi and /i-, for disk and bulge, and 
the fraction of the total light coming from the bulge. 





Disk 




Bulge 




Lb/(Lb + id) 




hi (kpc) 


h, (kpc) 


hi (kpc) 


h, (kpc) 




B 


5.8 


0.3 


1.2 


0.75 


0.39 


I 


5.0 


0.7 


1.2 


0.75 


0.45 



where Cbi = I miAB) 1 1 miAI) is the intrinsic I - B colour of the 
stellar population. An improved estimate of Jinti (/) is obtained 
by transforming the extinction factor from B to I via the stan- 
dard reddening law (ISavage&Mathislll979l) . The calculation 
is repeated several times 

We approximate xb and xi from a simple model of the 
galaxy wherein the intrinsic (i.e., three-dimensional) light dis- 
tribution 77 (x) for the disk and bulge are taken to be double 
exponentials. In cylindrical coordinates for M31, we have 



//' (x) = 77oe-'/"'<e 



-rlli,„-zlh[ 



(17) 



where the superscript / denotes either the disk or bulge, tjq 
is a normalisation constant, and and h, are the radial and 
vertical scale lengths, respectively. Different scale lengths are 
used for B and 1 because the two bands have different sen- 
sitivities to young and old populations of stars. Young stars 
tend to lie closer to the disk mid-plane than old ones. Our 
choices for the parameters are given in Table |3 The val- 
ues of the disk scale lengths and the bulg e-to-disk-ratios are 
taken from IWalterbos & KennicuttI (Il988l) The scale lengths 
for bulge are adapted from their de Vaucouleurs fit while the 
disk scale heights are based on the distribution of different stel- 
lar populations in the Milky Way disk. Th e obser vables Jobs(^) 
and Jobs(fi) are from 'Guhathakurt a et alJ ( l2005l) who cover a 
1.7°x5° field centred on M31. We derive colour profiles from 
th eir mosaics which are found to be similar to the profiles 
in IWal terbos & Kennicutt (1988). The colour is approximately 
constant within 30" and becomes bluer at larger radii. 

Our I-band extinction map for M3 1 is shown in Fig.^] The 
major dust lanes are clearly visible in the northern field and, as 
expected, the derived extinction is much larger on the near side 
of the galaxy than on the far side. The 1-band attenuation is 
< 40% and reaches a maximum in the innermost dust lane and 
a few smaller complexes. 
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Our model almost certainly underestimates the effect of 
extinction across the M31 disk. The approximation Jobs(^) - 
-^intr(^) is a poor starting point in the limit of large optical 
depths. For t » 1, most of the light in both B and 1 from be- 
hind the dust layer is absorbed and therefore IobsiB)/IobsiI) - 
Cbi- However substituting this result into equation 1151 gives 
exp(-T) ^ 1, an obvious contradiction. By the same token, 
if the dust is distributed in high-r clumps, then / and B wave- 
lengths will be absorbed by equal amounts given essentially 
by the geometric cross section of the clumps. Moreover, the 
thin-layer approxi mation tends to yield an under estimate of the 
extinction factor dWalterbos & Kennicuttlfl988i) . Finally, scat- 
tering increases the flux observed towards the dust lanes and 
therefore also leads one to underestimate the extinction factor 
Some of these problems can be solved by using infrared data 
in the construction of the extinction map. In a future paper we 
plan to use 2MASS data in order to derive a more accurate 
model for differential extinction in M3 1 . 

We can use the distribution of variable stars in our survey 
to test and refine the extinction model. The underlying assump- 



tion of this exercise is that the intrinsic distribution of vari- 
ables is the same on the near and far sides of the disk. We 
begin by determining the perio ds of the variable stars using 
a multi-harmonic periodogram ( Schwarzenberg- Czernvll996h 
suitably modified to allow for unevenly sampled data. A six- 
term Fourier series is then fit to each lightcurve yielding addi- 
tional information such as the amplitude of the flux variations. 
Only variables with lightcurves that are well-fit by the Fourier 
series are used. 

We will use LPVs to test the extinction model because 
they generally belong to quite old stellar populations. This is 
an advantage because the majority of the microlensing source 
stars also belong to older populations which are more smoothly 
distributed over the galaxy than younger variables such as 
Cepheids. We select LPVs with periods between 150 and 650 
days and focus on two regions of our INT fields. One of these is 
located on the near-side of the disk where extinction is expected 
to be high while the other is located symmetrically about the 
M31 centre on the far side. Fig.[21shows the spatial distribu- 
tion of the LPVs. Since extinction reduces the amplitude of the 
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Fig. 15. The distribution of the LPVs in M3 1 with the two sym- 
metrically placed regions used for the LPV amplitude analysis 
indicated. The northern field is located on the near side and 
contains some of the most heavily extincted parts, the southern 
field is on the far side and hardly affected by extinction. Thes e 
regions are similar to N2 and S2 regions from lAn et al.ll2004al) . 
only adjusted to avoid the part of the southern INT field that is 
not used in our analysis. 



flux variations and the average flux by the same factor we can 
study extinction by comparing the distributions in AF for the 
near and far sides. These flux variation distributions are shown 
in Fig. For low AF, where the shapes of the distributions 
are dominated by the detection efficiency, results for the near 
and far side agree. For high AF, where the detection efliciency 
for variables approaches 100%, one finds a large discrepancy 
between the near and far-side distributions. 

To test whether this discrepancy is indeed due to extinction 
we transform the coordinates of LPVs on the far side to their 
mirror image on the near side. The amplitude of the flux vari- 
ation is then reduced by the model extinction factor suitably 
transformed from 1 to r' ( Sav age & Mathis 1979). The new dis- 
tribution, shown in Fig.^] is still significantly above the near- 
side distribution at large AF though it does provide a better 
match than the original far-side distribution. The implication is 
that our model underestimates extinction. To explore this point 
further we consider models in which t is replaced by cr where 
c > l.InFig.lH we show the distributions of the far side LPVs 
for T — > 2r (long-dashed line) and r — » 2.5t (dot-dashed line). 
Apparently, the bright end of the (mirror) far-side distribution 
with T increased by a factor of 2.5 agrees with the bright end of 
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Fig. 16. Luminosity functions of LPVs in the 2 symmetri- 
cally placed regions. The far side flux distributions were scaled 
slightly to correct for small differences in area due to the gaps 
between the CCDs. The solid line is for the near side region and 
the dotted for the uncorrected far side region. The short-dashed, 
long-dashed, and dot-dashed lines are far side distributions cor- 
rected for increasing levels of extinction. 

the near-side distribution. We therefore conclude that our orig- 
inal model does indeed underestimate the effects of extinction. 
In some places this will be stronger than in others, but over the 
probed region the model underestimates extinction effectively 
by perhaps a factor of 2.5 in r. 

8. Theoretical predictions 

The detection efficiencies found in Sect. |6l allow us to predict 
the number and distribution of events given a specific model for 
the galaxy. Though M3 1 is one of the best studied galaxies, a 
number of the parameters crucial for microlensing calculations, 
are not well-known. Chief among these are the mass-to-light 
ratios of the disk and bulge, (M/L)^ and (MIL\, respectively. 
The light distributions for these components are constrained by 
the surface brightness profile while the mass distributions of the 
disk, bulge, and halo are constrained by the rotation curve and 
line-of-sight velocity dispersion profile. However, the mass-to- 
light ratios are poorly constrained primarily because the shapes 
of the disk and halo contributions to the rotation curve are simi- 
lar (e.g. jan Albada et al. 1985). One can compensate for an in- 
crease in (M/L)(j by dec reasing the overall den sity of the halo. 
Stellar synthesis models teell & de ■Tongl200 1^■ combined with 
observations of the colour profile of M31, can be used to con- 
strain the mass-to-light ratios though these models come with 
their own internal scatter and assumptions. Another poorly con- 
strained parameter is the thickness of the disk which affects the 
disk-disk self-lensing rate. 

In this section we describe theoretical calculations for the 
expected number of events in the MEGA-INT survey. We con- 
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sider a suite of M3 1 models which span a wide range of values 
in (M/L)(j and {M/L\. The dependence of the microlensing 
rate on other parameters is also explored. 

8.1. Self-consistent models of M31 

The standard practice for modelling disk galaxies is to choose 
simple functional forms for the space density of the disk, bulge, 
and halo tuned to fit observational data. For microlensing calcu- 
lations, velocity distributions are also required. Typically, one 
assumes that the velocity distribution for each of the compo- 
nents is isotropic, isothermal, and Maxwellian with a disper- 
sion given by the depth of the gravitational potential or, in the 
case of the b ulge, the observed line-of-sight velocity disper- 
sion. (But see lKerins et alJ ( l200lb where the effects of velocity 
anisotropy are discussed.) This approach can lead to a variety 
of problems. First, these "mass models" do not necessarily rep- 
resent equilibrium configurations, that is, self-consistent solu- 
tions to the collisionless Boltzmann and Poisson equations. A 
system initially specified by the model may well relax to a very 
different state. Another issue concerns dynamical instability. 
Self-gravitating rotationally supported disks form strong bars. 
This instability may be weaker or absent altogether if the disk is 
supported, at least in part, by the bulge and/or halo. Therefore, 
models with very high (M/L)^ are the most susceptible to bar 
formation and can be ruled out. 

In order to overcome these difficulties we use new, 
multi-component models for disk galaxies developed by 
[widrow & Dubinski (2005). The models assume axisymmetry 
and inco rporate an exponential disk, a Hernquist mode l bulge 
jHernquisl .1990) . and an NFW halo dNavarro et alJ ^996.^ . 
They represent self-consistent equilibrium solutions to the cou- 
pled Poisson and collisionless Boltzman n equations and are 
genera ted using the approach described in lKuiiken & Dubinski 

The phase-space distribution functions (DFs) for the disk, 
bulge, and halo (/disk, /bulge, and /haio respectively) are chosen 
analytic functions of the integrals of motion. For the axisym- 
metric and time-independent system considered here, the an- 
gular momentum about the symmetry axis, J,, and the energy, 
E, are integrals of motion. Widrow & Dubinski (2005) assume 
that /halo depends only on the energy while /bulge incorporates 
a /j. -dependence into the Hernquist model DF to allow for ro- 
tation. For both h alo and bu lge, the DFs are "lowered" as with 
the King model ( lKingll 966) so that the density goes to zero at 
a finite "truncation" radius. The disk DF is a function of E, 7 , 
and an approximate third integral of motion, E^, which corre- 
sponds to t he energy associ ated with vertical motions of stars 
in the disk d Kuiiken & Dubinski. 1995^) . 

Self-consistency requires that the space density, p, and 
gravitational potential, tfr, satisfy the following two equations: 



and 



VV = 47rGp . 



(18) 



(19) 



Self-consistency is achieved through an iterative scheme and 
spherical harmonic expansion of p and tfr. Straightforward tech- 
niques allow one to generate an N-body representation suitable 
for pseudo-observations of the type described below. The N- 
body representations also provide very clean initial conditions 
for numerical simulations of bar formation and disk warping 
and heating. 

The DFs are described by 15 parameters which can 
be tuned to fit a wide range of observations. In addition, 
one must specify mass-to-light ratios if photometric data 
is used. Our strategy is to compare pseudo-observations of 
M31 with actual observational data to yield a ;t'^-statistic. 
Minimisation of over the model parameter space - per- 
formed in Widr ow & Dubin skil (l2005l) by the downhill simplex 
method (see e.g. Press et al. 1992) - lead s to a b est-fit model. 

Following .Widrow & Dubinskil ( l2005l) (see, also 
IWidrow et alJ ll2003[) who carried ou t a sin ular exercise 
with the original feuiiken & Dubinskil [[$95') models) we 
utilise measurements of the surface brightness profile, 
rotation curve, and inner (that is, bulge region) velocity 
profiles. We use R-band surface brig htness profiles for th e 
major and minor axe s from Walterbos & Kennicutll (Il988l) . 
(Widrow & Dubins kil tOO^ used the global su rface bright- 
ness profile from Walterbos & Kennicutt ( 1988) which was 
obtained by averaging the light distribution in elliptical rings. 
The use here of both major and minor axis profiles should yield 
a more faithful bulge-disk decomposition.) The theoretical 
profiles are corrected for internal extinction using the model 
described in the previous section. In addition, a correction for 
Galactic extinction is included. We assume photometric errors 
of 0.2 mag. We use a composite r otatio n curve constr ucted 
from observations by iKentI (Il989l) and iBraunI ( Il99ll) that 
run from 2 to 25 kpc in galactocentric radius. Values and 
error bars for the circular speed are obtained at intervals of 
lOarc min ^ 2.2 kpc using kernal smoothing dWidrow et alJ 
'2003'). Finally, we use kinematic measurements from McE lrovi 
( 1 983) to constrain the dynamics in the innermost part of the 
galaxy. We smooth his data along the minor axis to give values 
for the line-of-sight stellar rotation and velocity dispersion at 
0.5 kpc and 1 .0 kpc. The values at these radii are insensitive 
to the effects of a central supermassive object and reflect the 
dynamics of the bulge stars with little disk contamination 
(McEkoy 1983). An overall for the model is calculated by 
combining results from the three types of data. Photometric 
and kinematic data are given equal weight; the circular rotation 
curve measurements are weighted more heavily than the bulge 
velocity and dispersion measurements. To be precise, we use 

1 , 2 



X 



(20) 



where xl\,p^ /t'Lige' ^^'^ ^rc ^'"^ individual ;t'^-statistics for 
the photometric, bulge kinematics, and rotation curve measure- 
ments. 

Our reference model (model Al) is constructed with 
(M/L)d - 2.4 and (M/L\ = 3.6. These value s are motivated 
by the stellar population synthesis models of iBell & de .Tong l 
d200ll) . Along the far side of the minor axis, where the sur- 
face brightness profile is relatively free of extinction, the B- R 
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Fig. 17. Comparison of pseudo-observations of model Al to 
real observations. Upper panel: model surface brightness pro- 
files (solid line s) along the major and minor a xis compared to 
observations by Walterbos & Kennicutt ( 1988) (dots). For clar- 
ity the profiles are shifted down in steps of 2 magnitudes. From 
the top down the profiles correspond to: SW major axis, NE 
major axis, SE minor axis (far side), and NW minor axis (near 
side). Lower panel: mode l rotat i on cur ve (s olid line) and com- 
bined rotation curve from lKentI (Il989l) and lBraunI ( Il99lh . The 
three lower lines correspond to the contributions to the rotation 
curve of the bulge (dotted), disk (long dash) and halo (short 
dash). 

colour is 1.8 in the bu l ge r egion and 1.6 in the disk region 
IWalterbos & K ennicutt' (1988). A correction for Galactic ex- 
tinction brings these numbers down by 0.18. Substituting into 
the appropriate formula from Table 1 o f lBell&de.Tongl(l200ll) 
yield the mass-to-light ratios chosen for this model. In Fig. 
El we compare predictions for model Al with observations. 
Shown are the surface brightness profiles along major and mi- 
nor axes and the circular rotation curve. Not shown is the excel- 
lent agreement between model and observations for the stellar 
rotation and dispersion measurements in the bulge region. The 
reduced statistic for this model is 1.06 (see Table|5}- 

In model Al, the scale height of the disk was fixed to 
a value of 1.0 kpc. Note that our model uses a sech^-law 
for the vertical structure of the disk. A sech^-scale height of 
1 kpc is roughly equivalent to an exponential scale height of 
0.5 - 0.7 kpc. The observations used in this study do not pro- 
vide a tight constraint on the scale height of the di sk and so we 
appeal to observations of edge-on disk galaxies. iKregel et"d] 



(I2OO2I) studied correlations between the (exponential) vertical 
scale height and other structural parameters such as the radial 
scale height and asymptotic circular speed in a sample of 34 
edge-on spirals. Using these correlations we arrive at an ex- 
ponential scale height for M31 of 0.6 kpc with a fairly large 
scatter 

We also fix the disk truncation radius for this model to 
28 kpc w hich is at the high end of the range favoured in 
Ixregel et alJ ((20021) . Lower values appear to be inconsistent 
with the measured surface brightness profile. The remaining 
parameters for the disk, bulge, and halo DFs are varied in order 
to minimise 

Table |5] outlines other models considered in this paper. 
Models Bl-El explore the {MIL\ - (M/L)^ plane. The for 
these models are generally quite low, a reflection of the model 
degeneracy mentioned above. In these models, disk and bulge 
"mass" are traded off agains t halo mass. Previous investigations 
dWidrow & Dubinskill2005l) suggest that model El is unstable 
to the formation of a strong bar while the other models are sta- 
ble against bar formation or perhaps allow for a weak bar 

The aforementioned models used values for the extinction 
factor derived in Sect.0 As discussed in that section, there are 
a number of reasons to expect that this model underestimates 
the amount of extinction in M31. Indeed, our analysis of the 
near-far asymmetry in LPVs favours a higher optical depth by 
a factor of 2.5, that is, the substitution e"'" — > e"^-^'". For this 
reason, we consider a parallel sequence of models, A2-F2, with 
high extinction. Note that the for these models are as good 
as if not better than those for the corresponding low-extinction 
models. 

8.2. Event rate calculation 

The event rate is calculated by performing integrals over the 
lens and source distribution functions. The rate for lenses to 
enter the lensing tube of a single source is 

d^R = ^^^-p^2REV^dl,d\id/3 (21) 
Ml 

where /; is the DF for the lens population, // is the observer-lens 
distance (Dql in the notation of equation|2}, is the transverse 
velocity of the lens with respect to the observer-source line-of- 
sight, and Aii is the mass of the lens. In writing this equation, 
we assume all lenses have the same mass. 

For a distribution of sources described by the DF equa- 
tionl^is replaced by the following expression for the rate per 
unit solid angle 

dR rf,(l,,y,) /,(/., V,) 

= I 2Rfv \ 

dQ. J Ml (M/L),L, 

X dlid\i l] dhd\s dp (22) 

where Is is the observer-source distance, (M/L)s is the mass- 
to-light ratio of the source and is the source luminosity. (For 
the moment, we treat all sources as being identical.) 



Jelte T.A. de Jong et al.: MACHOs in M31? Absence of evidence but not evidence of absence 



17 



We perform the integrals using a Monte Carlo method. The 
DFs are sampled at discrete points: 



p i=i 



where p € {I, s), "Lp is the surface density of either lens 
or source population, and Np is the number of points used 
to Monte Carlo either lens or source populations. The nine- 
dimensional integral in equationl22lis replaced by a double sum 
and an integral over yS: 



where 

Ssl - 

and 



NiMiNsLs {MIL), 



(24) 



(25) 



(26) 



Note that S depends on the line of sight densities of the lens 
and source distributions along with characteristics of the two 
populations. depends on the coordinates and velocities of 
the lens and source (hence the ij subscripts). The sum is re- 
stricted to lens-source pairs with li < Is- For each lens-source 
pair, the Einstein crossing time, t^^ij is easily calculated. The 
differential event rate is then 



(fR 
dOdtE 



= Ssl J] J^" d/mij6(Mij - fE) . 



8.3. Stellar and MACHO populations 



The formulae in the previous section apply to the six lens- 
source combinations in our model: disk-disk, disk-bulge, 
bulge-disk, bulge-bulge, halo-disk, and halo-bulge. As written 
the formulae assume homogeneous populations. For the disk 
and bulge populations, we modify equation l27l to include in- 
tegrals over the mass and luminosity functions as appropriate. 
We write the luminosity function (LF) as 



dN 



dMR 

and the mass function as 
dN 

— = Bh{M,MQ) 
dM 



(28) 



(29) 



where A and B are normalisation constants and Mq is the 
lower bound for the mass funct i on (MF ). We take the func- 
tion g from iMamon & Soneiral (Il982h and the function h 
from lBinnev & MerrifieldVlQQSS (their equation 5.16) with the 
power-law form dNjdM oc Al"'^ extended to Mq. A and B are 
evaluated separately for the disk and bulge populations. In the 
case of the disk, we assume that 30% of the mass is in the form 
of gas. The LF is normalised to give L - Lq with the proviso 



that Ls in equation|52is given in solar units. To determine the 
normalisation constant B of the mass function, we write 



(23) Bh{Me,Ma} 



dN_dMv\ 
dMv dM j 



(30) 



where the V-band LF is again fro mlMamon & Soneiral ( Il982l) 
and dMv/dM is from Ki'oupa et'al.l (Il993l) . Equation BUI is 
evaluated at solar values for convenience. The relation 



M\ _ jBh(M,MQ)MdM 
J Ag{MR)L{MR)dMR 



(31) 



can then be solved for Mo- Thus, a disk with high M/ L contains 
more low-mass stars than a disk with low M/L. 

For simplicity, and because we lack a model for what 
MACHOs actually are, we assume all MACHOs have the same 
mass. Mm that is 



^^6{M-Mm) ■ 
dM 



(32) 



The value of Mm will directly determine the number den- 
sity of MACHOs for a given halo mass density. Since the 
MACHOs only provide lenses and no sources for microlens- 
ing, a higher value of Mm and thus a lower number density, 
will result in a lower number of microlensing events. A given 
value of Mm can be considered as the average mass of a more 
elaborate MACHO mass function. 



•^^^^ 8.4. Theoretical prediction for the number of events 



Recall that the efficiency e is written as a function of Jfwhm 
and AFmax- (The efficiency also depends on the Une of sight.) 
These quantities are explicit functions of /3, F^, and t^. Thus, 
the expected number of events per unit soUd angle is 



dMRg (Mr) 



X dMih(M,Mo)'Ru£(tFWHM, AF) 



(33) 



where E is the overall duration of the experiment. Our survey 
covers four half-year seasons and so, with our choice of units 
for e and dR/dO., we have E -2. 

The number of events expected in each of the 250 bins used 
for the extinction calculation and labelled by "k" is 



, dE\ 



(34) 



where tsD. - 9 arcmin is the angular area of a bin. Sk car- 
ries an additional label (suppressed for notational simplicity) 
which denotes the lens-source combination. The total number 
of events is & - Yj&k- 
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8.5. Binary lenses 

Our microlensing selection criteria are based on the assump- 
tion that the lenses are single point-mass objects. However, 
at least half of all stars are members of multiple star sys- 
tems. Microlensing lightcurves for a lens composed of two or 
more point masses can deviate significantly from the standard 
lightcurve (Schneider & Weiss 1986) and may therefore escape 
detection. The deviations are strongest when the source crosses 
or comes close to the so-called caustics, positions in the source 
plane where the magnification factor is formally infinite. (The 
actual magnification factor is finite due to the finite size of 
the source.) The size of the caustic region is largest when the 
separation of the components of the lens is comparable to the 
Einstein radius corresp onding to the total mass (equation |2j- 
iMao & Paczvnskj| (!l991) estimated that ~10% of microlensing 
events towards the bulge of the Milky Way (mainly self-lensing 
events) should show strong binary characteristics such as caus- 
tic crossings. Since the Einstein radius for bulge-bulge self- 
lensing toward the MiUcy Way and M31 ar e comparable, we 
can e xpect a similar 10% effect in our survev. lBaltz & Gondold 
1I2OO I) perform a similar analysis for pixel-lensing surveys and 
estimate that in the order of 6% of self-lensing events from 
normal stellar populations will exhibit caustic crossings. Since 
the majority of detected events will have low signal-to-noise, 
we can assume that deviations other than caustic crossings in 
most cases will not strongly affect our detection efficiency. 
Therefore, to account for binary lenses, the calculated theoret- 
ical predictions for self-lensing are revised downward by 10%. 

8.6. Results 

Table 5 presents the theoretical predictions for the total number 
of events expected in the MEGA-INT four-year survey. The 
results are given for both self-lensing (Sseif) and halo lensing 
(£haio)- The values quoted for fihaio assume 100% of the halo is 
in the form of MACHOs. In other words, these values should 
be multiplied by the MACHO halo fraction in order to get the 
expected number of events for a MACHO component. We note 
that lensing by the Milky Way halo is not included in these 
results. This possible contribution is expected to be small, since 
the number of microlensing events from a 100% MW halo is 
a few' times lower than fo r a 100% M31 halo (Gvuk & Crotti 
l2000HBaillon et alJl993b for MACHO masses around O.SMq. 

We also consider the near-far asymmetry for self and halo 
lensing. In Fig. ^] we show the cumulative distribution of 
events for self and halo lensing as a function of the distance 
from the major axis, s. We take s to be positive on the far side 
of the disk. For this plot, we choose model Al but since the dis- 
tributions are normalised to give 14 total events, the difference 
between the models is rather inconsequential. We see that both 
self and halo lensing models do a good job of describing the 
event distribution in the inner 0.2°. The halo distribution does 
a somewhat better job of modelling the three events between 
s - 0.2° and s - 0.3°. Neither halo nor self lensing models 
predict anywhere near two events for s > 0.35°. 



To further explore the distribution, we define the asymme- 
try parameter 

o 

In Table 5 we give values for J?lseif and ^haio- We also pro- 
vide an average yiave which assumes that MACHOs make up 
the shortfall between the expected number of events and the 
observed value of 14. In cases where the expected number of 
events is greater than 14, we set ^ave - -^seif ■ The asymmetry 
parameter for the 14 candidate events is ^data = 0.125. 

The general trend, in terms of total expected number of 
events, is that as the mass-to-light ratios are increased, £seif 
increases and £haio decreases. There are counter examples. In 
model CI, the (M/L\ (as compared with model Al) leads to a 
less massive disk and lower £seif ■ Recall that for each choice of 
mass-to-light ratios, the remaining parameters are adjusted to 
minimise The process can lead to rather complicated inter- 
dependencies between the model parameters. The self-lensing 
rate decreases with decreasing as illustrated with model Fl. 
The self-lensing rate is generally reduced in the high extinction 
models relative to the low extinction ones. Finally we see that 
the halo event rate decreases with increasing MACHO mass. 
Models G and H illustrate this point and span the range in Mm 
identified by Alcock et al. (2000.) as the most probable mass 
range for Milky Way MACHOs. 

The timescale distribution is easily calculated using the 
method outlined in the previous section. Essentially, one calcu- 
lates fpwHM for each lens-source pair in the Monte Carlo sum. 
In Fig.^]we show the cumulative timescale distribution of our 
candidate microlensing event sample and model Al. In con- 
structing the curves for self and halo lensing, we have scaled 
the distributions to give a total of 14 events. 

9. Discussion 

The numbers expected for events due to self-lensing across the 
models probed in Table |5] fall within the narrow range of 10- 
16. The relative insensitivity of fiseif to changes in the mass-to- 
light ratios is a result of our approach to constructing models; 
changes in {M/L\ and (M/L)^ are compensated by changes in 
the structural parameters of the disk, bulge, and halo so as to 
minimise for the fit to the rotation curve and surface bright- 
ness data. Consider models D 1 and E 1 . The mass-to-light ratios 
differ by a factor of ~ 2 while fiseif differs by only a factor of 
1.4; with the low M/L values in model Dl, the rotation curve 
data drive up the disk and bulge luminosity distributions at the 
expense of a poorer fit for the photometric data. A balance is 
struck and the net result is that the change in £seif is signifi- 
cantly smaller than what one might expect. 

The consistency of the number of candidate events with the 
number of predicted self-lensing events is contrary to the re- 
sults of the analysis of the first th ree seasons of INT data by 
the POINT- AGAPE collaboration, talchi Novati et al.' ("2005) 
present six high quality, short duration microlensing candidates 
with one of these events attributed to M32-M31 lensing. They 
also model the detection efficiency and calculate number of ex- 
pected self- and halo-lensing events for a variety of M3 1 mod- 
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Table 5. Results of the microlensing modelling using self-consistent M31 models. In the first columns some model parameters 
and the combined are listed. The remaining columns contain the predicted number of events due to self-lensing (Sseif), due to 
halo-lensing {&haio), the asymmetry of the self-lensing itself), of the halo-lensing {Jihaio), and of the combination of both {Jiave)- 
The number of self-lensing events &seif has been corrected for the fact that ~10% of the events will show strong binary effects 
and therefore be selected against. The microlensing event rate due to the halo &haio is for a 100% MACHO halo, i.e. all of the 
halo mass is assumed to be in the MACHOs. For calculating the combined self- and halo-lensing asymmetry parameter Jlave a 
smaller fraction of the halo mass is assumed to be in MACHOs, namely the amount necessary to make up the difference, if any, 
between &seif and the observed number of 14 candidate events. The disk scale heights h, are sech^ scale heights. The upper, low 
extinction part of the table contains models with internal extinction values as derived in section^ while the lower, high extinction 
part contains models with increased extinction, as motivated by our analysis of the LPV amplitudes. 

Low extinction 



Models with mMACHO=0.5 Mq and /i;=1.0 kpc 





(M/L)d 


(M/L\ 




fisclf 


fihalo 


'fisclf 


■fihalo 




Al 


2.4 


3.6 


1.06 


14.2 


30.9 


0.037 


0.086 


0.037 


Bl 


2.4 


2.9 


1.17 


13.4 


31.5 


0.031 


0.085 


0.033 


CI 


2.4 


4.3 


1.02 


13.1 


29.6 


0.039 


0.092 


0.043 


Dl 


1.8 


2.4 


1.34 


11.3 


35.5 


0.031 


0.082 


0.041 


El 


3.6 


4.4 


1.03 


15.8 


24.6 


0.030 


0.091 


0.030 


Models with (M/L)d=2.4 and (M/L)b=3.6 




hz 


Mm 


X- 


fisclf 


fihalo 


-^sclf 


■fihalo 


'-^ave 


Fl 


0.5 


0.5 


1.10 


12.5 


30.7 


0.037 


0.084 


0.042 


Gl 


1.0 


0.1 


1.06 


14.2 


43.1 


0.037 


0.088 


0.037 


HI 


1.0 


1.0 


1.06 


14.2 


25.9 


0.037 


0.085 


0.037 



High extinction 



Models with Mm = O.SMq and /i,= 1.0 kpc 





(M/L)d 


{MIL\ 


x' 


fisdf 


fihalo 


-fisclf 


-fihalo 


i^ave 


A2 


2.4 


3.6 


0.99 


12.4 


28.6 


0.052 


0.095 


0.057 


B2 


2.4 


2.9 


1.08 


12.2 


32.6 


0.046 


0.094 


0.052 


C2 


2.4 


4.3 


0.99 


14.5 


29.6 


0.056 


0.098 


0.056 


D2 


1.8 


2.4 


1.23 


10.3 


34.5 


0.045 


0.095 


0.058 


E2 


3.6 


4.4 


1.04 


14.2 


22.8 


0.046 


0.105 


0.046 


Models with (M/L)d=2.4 and (M/L)b=3.6 




hz 


Mm 


X' 


fisclf 


fihalo 


-^sclf 


-fihalo 


i^ave 


F2 


0.5 


0.5 


1.06 


11.2 


30.5 


0.052 


0.095 


0.061 


G2 


1.0 


0.1 


0.99 


12.4 


39.1 


0.052 


0.098 


0.057 


H2 


1.0 


1.0 


0.99 


12.4 


23.8 


0.052 


0.093 


0.057 



els. In all of their models, the number of events for self-lensing 
is predicted to be less than ~ 1.5. Since this number is sig- 
nificantly less than the observed number, they conclude that 
some of the events are due to MACHOs and estimate that the 

M ACHO halo fraction is at le ast 20%. 

ICalchi Novati et alJ ( l2005t) use the m odel fromlKerins et alJ 
1I2OOI') which features a bulge following Kent an expo- 

nential sech^ disk and a spherical, nearly isothermal halo. They 
use the same structu ral parameters for the three components as 
iKerins et alJ ( I2OOI ') but take (M/LB)h = 3 and {M/L)d = 4. 
This model for the stellar mass distribution in M31 predicts 
an inner rotation curve that is significantly lower than the ob- 
served one, and so an extra 'dark bulge' component is required 
as well as the isothermal halo. £alchi Novati et al. (2005) do 
not consider microlensing by this dark bulge in their model, 
but instead attribute all surplus microlensing to the halo. In our 
model the stellar bulge is more massive, with M/L that is suffi- 



cient to reproduce the inner rotation curve, and there is no non- 
lensing dark bulge component. It appears to provide sufficient 
microlensing events to explain the observations. 

Furthermore, the choice of 0.3 kpc for the sech^ scale 
height is small by perhaps a factor of 3 if M3 1 is a typical spi - 
ral galaxy as represented in the survey by Krege l et alJ 
Thickening the disk increases the disk-disk self-lensing rate. 

For our models, the number of events due to self-lensing 
is consistent with the total number of events observed but not 
inconsistent with a significant MACHO fraction for the halo of 
M3 1 . We can make this statement more quantitative by treating 
halo events as a Poisson process with bac kground due to self- 
lensing and employing the approach of iFeldman & Cousins! 
( 1998). We let n be the number of observed events consisting 
of MACHO events with mean /fihalo, where / is the MACHO 
fraction, and a background due to self-lensing with known 
mean fiseif. For this analysis, we ignore the background due 
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0.2 
s(degrees) 



Fig. 18. Cumulative event distribution as a function of distance 
from the major axis (in degrees). Shown are the data (dots), 
self-lensing distribution (solid line), and halo-lensing distribu- 
tion (dotted line). Both self- and halo-lensing lines are scaled 
to give a total of 14 events. 




tpwHM (days) 

Fig. 19. Cumulative microlensing event distribution as a func- 
tion of timescale. The line and point-types are the same as in 
Fig.lHl 



to variables and background supernovae. The probability dis- 
tribution function is 



P in\f) = (/£hdo + Sself)" exp [- (/fihalo + fiself)] /«! 

To obtain confidence intervals for /: 



(36) 



1. Calculate P(n\f) for values of / e {0, 1) and sort from 
high to low. The maximum of P defines the most probable 



value of /. The values of P are normalised so that the sum 
of all sampled values of f is 1 . 
2. Accept values of / starting from the highest value of P un- 
til the sum of P exceeds the desired confidence level. The 
largest and smallest values of accepted / define the confi- 
dence interval. 

In Table |6l we provide most probable values of / and 95% 
confidence intervals for all of the models in Table|5] We provide 
these values both for the case of the full sample of 14 observed 
candidate events (n=14), as well as for the case of 11 observed 
events (n=l 1), for reasons discussed below. 

We next turn to the distribution of events across the M31 
disk as represented by the asymmetry parameters. From Table 
|5]we see that ^seif < '^hMo < ^i-na- The (weak) asymmetry in 
the self-lensing distribution is due to extinction. Note that the 
values are significantly below J?ldata even for the high extinction 
models. 

The asymmetry parameter for the halo is significantly 
higher than that for self-lensing events and close to, though still 
below, yidata- However, the asymmetry parameter for combina- 
tions of self and halo lensing are well below .^data- Evidently, 
the distribution of candidate events is diflicult to explain with 
any reasonable combination of self and halo lensing. 

The large asymmetry in the data is due, for the most part, 
to events 11, 13, and 14 (see Table 0. It is therefore worth 
considering alternative explanations for these events. As argued 
in Paulin-Henriksson et al. (2002), the lens for event 11 likely 
resides in M32 and since we have not included M32 in our 
model, this event should be removed from the analysis. Doing 
so leads to a modest reduction in ^Idata- 

Events 13 and 14 may be more diflicult to explain. For 
model Al, the predicted number of self-lensing events with s > 
s (event 18) is 0.005 while the predicted number of MACHO 
events in the same range in s is 0.14/. Thus, the probability of 
having two events either from self or halo lensing is exceed- 
ingly small, unless the halo fraction is very large. However, 
since some contamination by variable stars of our sample can 
not be excluded, one or both of these events may be a variable 
star. We note, for example, that event 13 has the lowest S/N in 
our sample. The probability of having one event for MACHO 
lensing with / - 0.20 is ~ 3%, small, but not vanishingly so. 

A closer inspection of the model is also warranted. 
Recall that our models assume axisymmetry whereas M31 
exhibits a variety of non-axisymmetric features such as disk 
warping. This poiiit is il lustrated in the isophotal map by 
iHodge & Kennicutll ( ll982l) . From the map, one finds that event 
13 lies on the B=24 (R=22.6) contour while model Al predicts 
R=23.5. Thus, the model may in fact underestimate the surface 
brightness of the disk by a factor of 2, and hence the disk-disk 
self-lensing rate by a factor of 4. (The reason for the discrep- 
ancy is not completely clear The contours on the far side do 
appear to be "boxier" than those predicted by the model.) 

It is interesting to note that events 13 and 14 are coinci- 
de nt with the l ocatio n of the giant stellar stream discovered 
by llbata et aTl (l200lh . This stream runs across the southern 
INT field, approximately perpendicular to the major axis and 
over M32. Indeed, M32 may be the progenitor of the stream 
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Table 6. Most probable value and 95% confiden ce limits for the 
MACHO halo fraction / from the Feldman & CousinsI lll998l) 
analysis, for the full sample and the case without candidate 
events 11, 13, and 14. 







14 events 




11 events 


model 


/best 


conf. interval 


/best 


conf. interval 


Al 


0. 


[0,0.28] 


0. 


[0.,0.21] 


Bl 


0.02 


[0,0.29] 


0. 


[0.,0.22] 


CI 


0.03 


[0,0.32] 


0. 


[0.,0.24] 


TM 


A AO 


rn n qai 
LU,U.30J 


A 

0. 


FA A TO! 
[V.,U.ZZ\ 


El 


0. 


[0,0.32] 


0. 


[0.,0.25] 


Fl 


0.05 


[0,0.32] 


0. 


[0.,0.24] 


Gl 


0. 


[0,0.20] 


0. 


[0.,0.15] 


HI 


0. 


[0,0.34] 


0. 


[0.,0.25] 


A2 


0.06 


[0.,0.35] 


0. 


[0.,0.25] 


B2 


0.06 


[0.,0.31] 


0. 


[0.,0.23] 


C2 


0. 


[0.,0.29] 


0. 


[0.,0.22] 


D2 


0.11 


[0.,0.33] 


0.02 


[0.,0.24] 


E2 


0. 


[0.,0.39] 


0. 


[0.,0.29] 


F2 


0.09 


[0.,0.35] 


0. 


[0.,0.26] 


G2 


0.04 


[0.,0.25] 


0. 


[0.,0.18] 


H2 


0.07 


[0.,0.42] 


0. 


[0.,0.31] 



Table 7. Observed number of events and the asymmetry of 
their spatial distribution, shown for the full sample of 14 events 
and for cases where the probable M32 event (11) and candi- 
date events 13 and 14 are ignored. The quoted errors are Icr 
errors, determined with the bootstrap method. Also shown is 
the asymmetry for the long-period variable stars (LPVs). 



Events used 


^data 


•^data 


Full sample 


14 


0.125 ± 0.046 


without 11 


13 


0.120 ±0.049 


without 13, 14 


12 


0.076 ± 0.034 


without 11, 13, 14 


11 


0.066 ± 0.034 


LPVs 


20,864 


0.071 ±0.001 



jMerrett et alJ2003l) . The average V-band sur face brightness of 
the stream is 30 + 0.5 mag arcsec"^ ('ibata e t aPbOOlh 

but this is measured far from the projected positions of events 
13 and 14. The surface brightness of the stream might be sig- 
nificantly higher near the position of M32. Perhaps the most 
conservative statement one can make about the stream is that it 
is not bright enough to distort the contours near events 13 and 
14, that is, it cannot be brighter than the disk at these radii. The 
microlensing event rate due to stars in the stream is of course 
enhanced relative to the rate for self-lensing by the ratio of the 
distance from the stream to the disk and the thickness of the 
disk, that is, by a factor of ~ 20. The stream-disk lensing rate 
might be further enhanced if the stars in the stream have a large 
proper motion relative to the disk. These arguments suggest 
that the number of stream-disk events in the vicinity of M32 
might be 0.03 -0.1; perhaps high enough to explain one event. 

Fig-EHlprovides a summary of our results with respect to 
the expected number of events and the asymmetry parameter 
The points with error bars represent the data for the 4 cases 
considered in Table The solid circles and lines correspond 
to the high extinction case; the open circles and dotted fines 
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Events 

Fig. 20. Asymmetries and event numbers for data and models. 
The points with error bars on the left show fi and for the 
sample of candidate events and the subsamples listed in table 
The solid lines correspond to the high and the dotted lines to 
the low extinction models from table |5l The dots are the pure 
self-lensing points, with the MACHO mass fraction increasing 
along the line. MACHO fractions of 20% are indicated with 
vertical lines. 

correspond to the low extinction case. The circles assume pure 
self lensing while the lines trace out the values for increasing 
MACHO fraction with the tick-mark indicating the position of 
/ - 0.2. Once again, we see that the asymmetry parameter for 
the data is higher than that for any of the models. Removing 
events 13 and 14 does improve the situation as does increasing 
the optical depth r; the asymmetry remains a little higher but 
consistent with the models. 

10. Conclusions 

This paper presents the analysis of four seasons of M3 1 ob- 
servations at the INT, a subset of the MEGA survey of M3 1 . 
The observations were carried out to search for MACHOs in 
the halo of M3 1 . Our fully automated search algorithm iden- 
tified 14 candidate microlensing events from over lO*" variable 
sources. Three of the candidates were previously unpublished. 
The spatial and timescale distributions are consistent with mi- 
crolensing. 

The core of this paper is the comparison of this candidate 
event sample with a calculation of the expected number of 
events from self and halo lensing. This calculation breaks into 
three parts: a model for the extinction across the M31 disk; a 
model for the detection efficiency; and a suite of self-consistent 
disk-bulge-halo models for M3 1 . 

The results with regard to the fundamental question of 
whether there is a significant MACHO fraction in the halo are 
inconclusive. Based on the total number of events, we find that 
the most probable MACHO halo fraction / varies between 
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and 0. 1 depending on the model. Our event rate analysis is con- 
sistent with a total absence of MACHOs as the confidence in- 
tervals for all of our models include f - 0- On the other hand 
we can not exclude some MACHO component, since the con- 
fidence intervals extend typically up to / = 0.25 and even up 
to / = 0.4 for a few models. 

The spatial distribution of the candidate events is highly 
asymmetric and does seem to favour a MACHO component. 
However, for different reasons it is questionable whether the 3 
candidate events that largely determine the asymmetry signal 
should be used in this analysis. Thus, we conclude that both 
from the observed number of events, and from their spatial dis- 
tribution we find no compelling evidence for the presence of 
MACHOs in the halo of M3 1 . 
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Appendix A: Candidate event lightcurves 

On the following pages, for each of the 14 candidate microlens- 
ing events in our sample, the r' and i' lightcurves and thumb- 
nails taken from the difference centred on the event positions 
are shown, together with a short discussion. Apart from the INT 
r' and i' data, KP4m R and I data points are also plotted in the 
lightcurves. The fits shown are however the fits done to only 
the INT data. 
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MEGA-ML-1 

Located close to the centre of M3 1 , this event has a rather noisy 
basehne. Apart from the background of very faint variables 
there are some variable sources clearly visible in the difference 
images. As can be seen in the thumbnails in figure lAJT b) a 
bright variable is located just a few pixels from the position of 
the candidate event. Another, fainter variable is seen at a simi- 
lar distance above and to the left. The other variable sources are 
further away and should have no influence on the photometry. 

MEGA-ML-2 

This candidate event is located very close to MEGA-ML-1 and 
therefore has the same problems connected to being close the 
centre of M31. In the thumbnails of days 94, 754, and 1208 
we see a variable source a few pixels to the left of the event 
position. This variable is brighter in r' than in i', which causes 
the r' baseline to be the most noisy. 

MEGA-ML-3 

This candidate event is also located close to the M31 centre. 
In figure |6l we already demonstrated that a very faint variable 
source is positioned ~0.25"away from this candidate event. In 
the i' thumbnails another variable is visible just above and to 
the right of the event. This variable has a bright episode be- 
tween days 440 and 480, causing the bump in the baseline in 
the i' lightcurve. 

MEGA-ML-7 

By far the brightest event in our sample, the thumbnails of 
MEGA-ML-7 show a very bright residual close to the peak 
centre. Since the peak occurs during the first season, some of 
the exposures used for creating the reference image contained 
a significant amount of the magnified flux, so that the base- 
line Ues at a negative diff'erence flux. There are some variables 
nearby, but none of them are close or bright enough to signif- 
icantly influence the photometry. The distance to the centre of 
M3 1 is also quite large (~22'), reducin g the background of faint 
variab le sources. As pointed out bv IPaulin-Henriksson et alJ 
( I2OO3I) . there are some systematic deviations from the best fit 
microlensing model. lAnetal'J l f20043) find that this anomaly 
can be explained by a binary lens. 

MEGA-ML-8 

This near side event is located ~23' from the centre of M31. 
A variable that is particularly bright in i' is situated about 
2.4" NW of the candidate event, but should not have much 
of an effect on the photometry. The baselines of the lightcurves 
indeed look stable and well-behaved. 

MEGA-ML-9 

Peak coverage is poor for this candidate event, but the baselines 
are stable. The thumbnails show quite a lot of faint variables, 



two of which are located very close, approximately 1" to the 
left of the event position, accounting for the noise in the i' base- 
line that is higher than in the r' lightcurve. 

MEGA-ML-1 

This event is a beautiful example of a combined lightcurve with 
KP4m and INT data. Peak coverage in INT i' is poor, but the 
KP4m I data points follow the fit (derived only from INT data) 
very well. A fairly bright variable is situated slightly above and 
to the right of the event position and there is a hint of a very 
faint variable about 1" to the left. Although the INT baseline 
in i' is noisy, the r' and both KP4m R and I lightcurves show an 
very stable and well-behaved baseline. 
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Fig. A.l. (a) Event 1: lightcurves. The two upper panels show the full r' and i' lightcurves of the microlensing event. In the lower 
left corner are zooms on the peak region. In the lower right corner the r' flux is plotted versus the i' flux; if the colour is constant, 
the points should lie on a straight line. Also drawn is the best fit microlensing model. The solid circles are points from the INT 
data, the open circles are from the KP4m data. The start of the INT survey, August 1st 1999, is used as the zeropoint for the 
timescale. 
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Fig. A.l. (b) Event 1: thumbnails. The two upper rows of thumbnails show are taken from r' and i' difference images during the 
peak of the candidate event. Selected thumbnails from the baseline are also shown in the two bottom rows. Each thumbnail is 
30x30 pixels or 10x10" in size. 
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Fig. A.2. (a) Event 2; lightcurves. See caption of lA.lf a). 
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Fig. A.2. (b) Event 2: thumbnails. See caption of lA.ir bt. 
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Fig. A.3. (a) Event 3: lightcurves. See caption of lA.lf a). 
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Fig. A.3. (b) Event 3: thumbnails. See caption of lA.lf bt. 
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Fig. A.4. (a) Event 7; lightcurves. See caption of lA.lf a). 



figapp84b.jpg 



Fig. A.4. (b) Event 7: thumbnails. See caption of lA.ll bt. 



figappQ5a.gif 



Fig. A.5. (a) Event 8: lightcurves. See caption of lA.lf a). 
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figapp05b.jpg 



Fig. A.5. (b) Event 8: thumbnails. See caption of lA.ir bl. 
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Fig. A.6. (a) Event 9: lightcurves. See caption of lA.lf a). 
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Fig. A.6. (b) Event 9; thumbnails. See caption of lA.lf b). 
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Fig. A.7. (a) Event 10: lightcurves. See caption of lA.lT a). 



figapp87b.jpg 



Fig. A.7. (b) Event 10: thumbnails. See caption of lA.lT b). 
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MEGA-ML-11 

A high signal-to-noise event with a good fit and stable base- 
line. There is some noise in the i' baseline, caused by the vari- 
able source that is visible in the thumbnails of days 6 and 756 
at ~1.3" above the event position. During the fourth observ- 
ing season a few bad columns were lying exactly on top of the 
event position, so that there is only 1 INT data point available. 
However, the KP4m data show that the baseUne remains flat 
everywhere. 

MEGA-ML-13 

This candidate event has the lowest signal-to-noise of our sam- 
ple. It is situated far out in the far side of the disk at ~3r from 
the centre of the galaxy and the relatively low galaxy back- 
ground makes it possible to detect these kind of faint events. 
Due to the y-axis scale the i' the baseline looks quite noisy, 
but it is in fact not significantly more so than for other candi- 
date events. The thumbnails of days 398 and 520 show that the 
closest variable source is located ~IA" below and to the left 
of the event, which explains the scatter in the i' baseUne. 



MEGA-ML-17 

The i' baseline is slightly noisy, but the r' and both KP4m 
lightcurves are well-behaved. In the thumbnails no very close 
variables are visible. 

MEGA-ML-18 

This candidate event shows quite large scatter in the baseline 
and also in the peak. Faint variables might be the culprits, al- 
though the event is not located very close to the galaxy centre 
(~15.1')- The thumbnails show no variable sources very close 
to the event position, however they do show that this event is 
situated on the edge of a fringe running diagonally across the 
thumbnails. This fringe and the fact that it can change posi- 
tion slightly between frames is the most probable cause for the 
noisy i' photometry. 



MEGA-ML-14 

At ~35.5' from the M3 1 centre, this candidate event is the most 
far out in the disk of all events in our sample. The i' photometry 
of this candidate event is compromised by the variable source 
at ~1.3". From the i' thumbnails one can also see that the event 
lies at the edge of a fringe, making the background in the lower 
half of the thumbnails brighter than in the upper half. This can 
also cause some extra scatter in the photometry. Overall, how- 
ever, the microlensing fit is very good and both INT and KP4m 
lightcurves show a stable baseline. 



MEGA-ML-15 

This event is again located close to the centre of M3 1 and pre- 
sumably has a strong background of faint variable sources. In 
the thumbnails also several variables are visible very close to 
the event position, both in r' and in i'. The lightcurve base- 
lines are rather noisy because of this, but show no coherent 
secondary bumps and the KP4m basehnes are very stable. 



MEGA-ML-16 

Not selec ted in our fir st analysis of the first two seasons of 
INT data (|de Jong et al.>.200 4) due to baseline variability, the i' 
lightcurve of this event is strongly influenced by a bright vari- 
able situated just 1.1" to the north. Using a smaller extraction 
aperture for the photometry in the present analysis, the i' base- 
Une is still very noisy and the same is true for the KP4m I-band 
data. The INT r' and KP4m R data are much better behaved 
and the r' peak is fit very well by the microlensing fit. 
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Fig. A.8. (a) Event 11: lightcurves. See caption of lA.lT a'). 
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Fig. A.8. (b) Event 11: thumbnails. See caption of lA.lf b). 
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Fig. A.9. (a) Event 13: lightcurves. See caption of lA.lf al. 



figapp09b.jpg 



Fig. A.9. (b) Event 13: thumbnails. See caption of lA.lf b). 
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Fig. A.IO. (a) Event 14: lightcurves. See caption of lA.lf at. 
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figapplQb. jpg 



Fig. A.IO. (b) Event 14: thumbnails. See caption of lA.ir b). 
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Fig. A.ll. (a) Event 15: lightcurves. See caption of lA.lf a). 
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Fig. A.ll. (b) Event 15: thumbnails. See caption of lA.lf bl. 
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Fig. A.12. (a) Event 16: lightcurves. See caption of lA.lf a'). 
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Fig. A.12. (b) Event 16: thumbnails. See caption of lA.lf bl. 
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Fig. A.13. (a) Event 17: lightcurves. See caption of lA.ir at. 



figappl3b.jpg 



Fig. A.13. (b) Event 17: thumbnails. See caption of lA.lf b). 
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Fig. A.14. (a) Event 18: lightcurves. See caption of lA.lf al. 



figappl8b.jpg 



Fig. A.14. (b) Event 18: thumbnails. See caption of lA.lf b). 
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